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We study the effect of quenched disorder on nonequiUbrium systems of interacting particles, 
specifically, driven diffusive lattice gases with spatially disordered jump rates. The exact steady- 
state measure is found for a class of models evolving by drop-push dynamics, allowing several physical 
quantities to be calculated. Dynamical correlations are studied numerically in one dimension. We 
conjecture that the relevance of quenched disorder depends crucially upon the speed of the kine- 
matic waves in the system. Time-dependent correlation functions, which monitor the dissipation 
of kinematic waves, behave as in pure system if the wave speed is non-zero. When the wave speed 
vanishes, e.g. for the disordered exclusion process close to half filling, disorder is strongly relevant 
and induces separation of phases with different macroscopic densities. In this case the exponent 
characterizing the dynamical correlation function changes. 
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What is the effect of quenched disorder on driven, 
nonequilibrium systems? This question is important in 
a number of physical situations involving flow in random 
media [Q . Theoretically, our understanding of these phe- 
nomena is based largely on numerical simulations and on 
the analysis of continuum equations for coarse-grained 
variables. In this Letter we obtain the exact steady state 
and static correlation function for a class of models of 
driven, interacting particles on a lattice - the drop-push 
process - with quenched site disorder in the hopping 
rates. Further, we study the time-dependence of hydro- 
dynamic fluctuations for this system in one dimension, 
and also for the disordered asymmetric exclusion process. 
We find that the behaviour of a system with current Jq 
and density p is largely determined by cq = dJo/dp. If 
the density is uniform on a macroscopic scale, cq is the 
mean speed of kinematic waves which transport density 
fluctuations through the system If Cq is nonzero, 
we conjecture that quenched disorder does not affect the 
asymptotic behaviour of the decay of fluctuations, and 
support this with extensive numerical results. By con- 
trast, we find that vanishing cq can signal the onset of 
disorder-induced phase separation with coexisting macro- 
scopic regions of different density, in which case the dy- 
namical behaviour is different. 

In a coarse-grained description of a l-d disordered 
current-carrying system, the nonuniform steady state 
density profile is described by a function po{x). The 
evolution of density fluctuations p = p{x,t) — po{x) is 
described phcnomcnologically by a stochastic evolution 
equation with spatially random coefficients: 

dtp = dx[D{x)dxp - c{x)p - \{x)p^... - ri{x,t)]. (1) 

This follows from the continuity equation dtp{x, t) + 
dxJ{x,t) = on writing the current as J{x,t) — 
Jsys(x, t) — D(x)dxP + ri{x, t), where D{x) is the space- 
dependent diffusion constant and rj is white noise; the 
systematic part of the current Jsys is expanded as Jq -I- 
c{x)£ + \{x)p^ .... The most relevant source of disorder 
in is the term with coefficient c{x) which represents 



the space-varying local speed of the kinematic wave of 
density fluctuations. The problem is equivalent to that 
of a moving interface in the presence of a certain type of 
columnar disorder; the interface height h{x,t) is related 
to p by /5 = dxh. Equation (1) then becomes 

dth = D{x)dxxh - c{x)dxh ~ \{x){dxhf... - ri{x,t). (2) 

The random coefficients c{x), X{x) now represent colum- 
nar disorder in the 2-d h-x space. It is important to 
understand how disorder affects scaling properties and 
to identify factors responsible for different universality 
classes ||] . A crucial difference between (|^) and the model 
studied by Krug Q is the absence of an additive quenched 
noise e{x) which models a frozen random contribution 
to interface mobility. Such a term, which strongly influ- 
ences static and dynamic properties, is necessarily absent 
in (|^) as a consequence of the spatial constancy of (Jsys) 
implied by particle conservation. Further, our model is 
distinct from depinning- threshold interface models with 
quenched point disorder . We show that the roughness 
exponent a — 1/2 in contrast to [|j and The decay 
of the kinematic wave in time leads to a power-law t^^ 
growth of the correlation function 

S'^{t) = {[h{x + Cat, t) - h{x, 0) - Jat]'^). (3) 

We conjecture that as long as cq is nonzero, leading power 
law behaviours are the same as in the pure system with no 
disorder (/? = i) This is supported by Monte Carlo 

simulation results for the disordered drop-push process, 
and for the disordered exclusion process in the regime 
|p - 1/2| > A with A 7^ 0. For \p-l/2\ < A we present 
evidence that Jq is independent of p. In this regime, 
Co vanishes and shocks separate macroscopic regions of 
different mean densities. 

The disordered drop-push process (DDPP) is a model 
of driven transport of carriers trapped in local regions of 
space, with each move involving a cascade of transfers 
through flUed traps (|,^. In l-d, on every site i of an 
L— site ring is a well which can accommodate at most U 
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Fig. 1 Typical DDPP configuration and move. 

particles (Fig. 1). Each well depth is chosen indepen- 
dently from a probability distribution P{li). The con- 
figuration of the system is specified by the set of parti- 
cle occupation numbers {^i} of all the wells. The dy- 
namics is stochastic, and in a small time interval dt, 
there is a probability e{ni\li)dt that a particle from well 
i hops out, and drops into well i + 1. If well i -I- 1 is 
already full, a particle from this well now gets pushed 
further right, and so on (Fig. 1). The cascade of 
adjacent-site jumps terminates once a particle drops into 
a well (say i -\- r) which was not completely full ear- 
lier. This elementary move thus changes configuration 

C = {.... Ui li^i ■■■h+r-l ni+r ■■■■} tO C = {.... 

Tii — 1 Zi+i ■■■U+r-1 TLi+r -\- 1 ....}. The jump ratcs 
e{ni\li) are pre-specified, and depend both on the depth 
of the well and the occupation. The rates and well depths 
are quenched random variables. 

The probability for the system to be in configuration C 
satisfies the standard master equation |^ with the tran- 
sition probabilities W[C — > C) identified with e's above. 
In the steady state, the total flux out of C equals the to- 
tal incoming flux. This is ensured if for every C obtained 
from C by an elementary transition there is a unique con- 
figuration C" such that in the steady state 



w{c -> c')^i{c) = w{c" ^ c)/i(c"). 



(4) 



Here /i(C) is the invariant measure and (|^) is the con- 
dition of pairwise balance We define weights for 
single-site occupations by Ui(0) = 1 and Uiijii) — 
n{l) n{2) ■ ■ ■ Ti{ni) if < ni < U. with Tijii) = 
eo/e("-i|'i), where eg is a microscopic rate. The measure 
for configuration C = {ui} has the product form 



(5) 



To show that this satisfies (^, we construct the config- 
uration C" corresponding to a given C and C as follows. 
Suppose the transition C — > C involves hopping a parti- 
cle at well i to well i + r with all wells in between full. 
Also suppose well i — s \s not full but all wells between 
i — s and i are full (Fig. 1). Configuration C" is identical 
to C at all sites except at the sites i — s and i, at which 
n'l_g = rii-s + I, n'l = rii — 1. Then (^ is satisfied, in 
view of the measure (pj). Since the dynamics is ergodic 
the invariant measure (ra) is unique for a fixed number of 
particles Q. This generahzes the result obtained earlier 
for the non-disordered case |q] , and is the first instance of 



an exact determination of the steady state of a spatially 
disordered nonequilibrium system of interacting particles 

The product measure form (|5|) allows us to calculate 
several physical quantities: the site densities, height- 
height correlation function, current Jqi and wave speed 
cq. In the limit of large L, it is convenient to in- 
troduce the generating function Zi = X^n =o 
where z is the fugacity. The probability of occurrence 
of configuration C is V{C) = tJi\c)z'^'' /]\^ Z,, where Np 
is the number of particles in the configuration. The 
steady state is characterized by uniform z, but inhomo- 
geneous site densities (n^) = zdlnZi/dz. The height- 
height correlation function r^(r) = {{hi+r — ^i)^) be- 
comes ^ ("-i)^) and can be evaluated. Dis- 



order averaging gives r(r) 
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, implying a = 1/2. 



To find the steady state current Jq, note that the cur- 
rent Ji,i+i across bond (i, i -f- 1) comes from jumps which 
either (a) originate from sites to the left of site i (with 
in-between wells full), or (h) originate from site i itself 
(a contribution ji). Class (a) events evidently also con- 
tribute to the current across bond — Of all events 
that contribute to Ji-is, class (a) is that subset of events 
in which site i is fully occupied. Since the probability of a 
jump between i — 1 and i is independent of the probability 
of occupation of site i, we have Ji^i+i = Ji-i,i Piih) + ji 
where ji = e('^i|^i) Pi{ni) = eo z (1 - Piik)) and 

Pi{ni) = Ui{ni)z^' /Zi. Noting that Ji_i,i = Ji,i+i = Jo, 
we find Jq — €q z. Since both Jq and p are known in 
terms of z, the macroscopic speed cq ~ dJ()/dp of the 
kinematic wave can be determined . 

Interestingly, the steady state measure and current can 
also be found in c? > 1, for models in which the ratio of 
hopping rates in different directions is independent of site 
and configuration, and the cascade of adjacent-site over- 
fiows in a single move is in the same direction | [Tl[ |. 

The other model we investigate is the disordered fully 
asymmetric simple exclusion process (DASEP) on a 1-d 
ring. In this model, each site can hold or 1 particle. A 
particle hop is attempted to the nearest-neighbour site 
on the right, and is completed only if the site in ques- 
tion is vacant. Attempt rates are associated with bonds 
and are disordered, with magnitudes varying from bond 
to bond, chosen from a binary distribution. Unlike in 
the DDPP, there is no analytical characterization of the 
steady state even with a single inhomogeneous bond. A 
numerical study has shown that if p is close to 1/2, a 
single weak bond acts as a blockage and produces un- 
equal densities over macroscopic length scales on either 
side of it, and a shock where the density profiles meet, far 
from the blockage ||l2|. We have studied the disordered 
case with an extensive number of weak bonds by Monte 
Carlo simulation and found that the steady state depends 
strongly on the filling. The current Jq varies smoothly 
with p provided \p — l/2\ > A (Regime A), whereas Jq 
has a single value if |p — 1/2| < A (Regime B) (Fig. 2a). 
The value of A depends on the ratio r of rates of the 
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Fig. 2 DASEP steady states (a) J vspfoir = f = 0.5 
(b) density profile for p — 0.80 (c) blowup of box in (b) 
(d) profile for p = 0.5 (e) blowup of box in (d). 

weak and strong bonds and on the fraction / of weak 
bonds, and is ~ 0.16 for r = / = 0.5. In Regime A, the 
density profile consists of a large number of shocks with 
a mean inter-shock spacing of a few lattice spacings (Fig. 
2c); the number of shocks scales with the system size. 
On length scales large compared to the inter-shock spac- 
ing, the density is roughly uniform. A semi-quantitative 
picture of the steady state can be obtained using a mean- 
field approximation, writing the current between sites i 
and i -|- 1 as J = ei^i+ipo{i){l — po{i + 1)) where the rate 
ei^i+i is eo (eo/2) for a strong (weak) bond. As J is the 
same in every bond, the densities {po{i)} satisfy a set 
of coupled, nonlinear equations, which can be iterated 
to convergence. The result is shown in Fig. 2c. The 
mean-field approximation evidently obtains locations of 
shocks fairly well, but not shapes of individual shocks. 
In Regime B, Monte Carlo results and mean-field calcu- 
lations show that the density is nonuniform on a macro- 
scopic scale (Fig. 2d), besides showing shock structure 
on the scale of a few lattice spacings (Fig. 2e). In 
this regime, the occurrence of long stretches of succes- 
sive weak bonds, coupled with the requirement of spatial 
uniformity of current, results in phase separation into 
high and low density phases, qualitatively as in the sin- 
gle defect case In both regimes, numerical results 
show r(r) ^ with a ~ 0.5 

Turning now to the dynamical behaviour of fluctua- 
tions in the steady state, we first consider systems which 
have a uniform density on macroscopic scales. The dy- 
namics of such 1-d systems is dominated by kinematic 
waves which transport density fluctuations through the 
system at a mean speed cq. Owing to quenched disorder. 



the local speed of the wave varies from one location to 
another (Fig. 3) and the question arises how wave dissi- 
pation, and thus S{t)j is affected. In the pure system, the 
evolution equation (|l|) has a:-independent coefficients and 
the long-time growth of S{t) is described by /? = 1/3 [0,1). 
In the disordered system, let us regard the coarse-grained 
medium as made of successive disordered patches, each 
independent of the other, and ask for the behaviour of 
a large-scale density fluctuation as it passes through the 
succession of patches. If Cq is nonzero, the probability of 
a density fluctuation in an infinite system returning to 
the same disorder patch dies down rapidly at long times, 
so it is a good approximation to regard the effect of dis- 
order as essentially uncorrelated in time. Further, since 
the speed c{x) is a spatially random function, the use of 
the averaged value cq in (j^) induces a noise of amplitude 
i^/^ in the location of the density fiuctuation. Since fluc- 
tuations in h scale as x^^^, the effect ('^ t^^^) on S"^ is 
subleading. Thus we conjecture that if cq 0, the long 
time behaviour of S'^{t) is ~ t^/'^, the same as for the pure 
system. Our argument differs from that used earlier for 
the effect of point disorder on unpinned, moving inter- 
faces p^ , as our case corresponds to columnar disorder 
in the interface language. The irrelevance of randomness 
in c{x) is consistent with straightforward power counting 
in (2). 



In our numerical determination of S{t) from simula- 
tions of the two types of lattice models, we defined h{i, t) 
as p{k,t) where io{t) is the position of a specific 

particle. We averaged results for S{t) over 40 time evo- 
lutions for a fixed realization of disorder, and then over 
several realizations. For drop-push dynamics, we con- 
sidered a model with only two types of wells A and B, 
distributed randomly. Each well can hold at most one 
particle, but the jump rates out of the two types of wells 
are different, say ea and cb- We used ca/^b = 0.5, and a 
fraction / = 0.5 of low rates. Since the placement of the 
wells is random, the essential feature of quenched disor- 
der is still present. Since (|) stiff holds, we start with 
an initial configuration of particles consistent with this 
product measure. The analysis is aided by the fact that 




Fig. 3 Typical time evolutions (a) DDPP (b) DASEP, 
p = 0.25 (c) DASEP, p = 0.5. Darker regions are particle- 
rich. The tilted streaks are kinematic waves. 
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Fig. 4 S^{t) for the DDPP and the DASEP for dif- 
ferent densities. Individual data sets have been shifted 
for clarity. The straight lines have slopes 0.84 and 0.67. 
Inset: The return time of the kinematic waves is given by 
the period of oscillations of the autocorrelation function. 



oppositely moving kinematic waves in different macro- 
scopic regions (Fig. 3c). 

We conclude with a recapitulation of our principal re- 
sults. For the DDPP, the steady state has a product 
measure form, and the current can be determined. For 
the DASEP, the steady state density profile shows many 
shocks and is quite well described by a mean-field approx- 
imation. Our conjecture, that disorder does not affect 
the dynamical universality class if the kinematic wave 
speed Co is nonzero, is borne out by simulation studies 
of the DDPP and the DASEP in Regime A. A vanishing 
Co can indicate phase coexistence with different densities 
in different macroscopic regions, as in Regime B of the 
DASEP; the dynamical behaviour is different in this case. 

We thank R. E. Amritkar, D. Dhar, S. Krishnamurthy 
and G. I. Mcnon for useful comments. 



Jo and Co are known explicitly. In the case of the DASEP, 
we chose r = / = 0.5 and used a system of size 8000, al- 
lowing it to relax for ~ 60, ODD Monte Carlo steps to 
achieve steady state. Further, Jo and co were estimated 
in two stages: First, a rough estimate of Jo was obtained 
by direct measurement, while Co (in Regime A) was es- 
timated from the return time of the kinematic wave in 
a smaller system (Fig. 4 inset). Then a more accurate 
estimate was obtained by minimizing S' in a large system 
with respect to Jo and co. Details of the minimization 
procedure will be given elsewhere . Fig. 4 shows data 
for two different densities for the DDPP and also Regime 
A of the DASEP, all corresponding to nonzero co. In 
all these cases, the data is consistent with ~ t'^^^, 
as for the pure system, supporting our conjecture. Fur- 
ther, we performed a direct check of the ~ t^^"^ growth 
of disorder-induced noise in locations of density fluctu- 
ations, discussed in the previous paragraph. Such noise 
should alter the growth law for the sliding-tag correla- 
tion function (which monitors fluctuations of tagged- 
particle locations) from t^/^^ to t^^"^ - a change we 
confirmed by simulation of the DDPP [ pl| . 

In Regime B of the DASEP, Sit) shows stronger fluctu- 
ations than in Regime A, from one realization of disorder 
to another. On averaging over 10 samples, we find S(€) 
grows as with /3 = 0.42 ± 0.02 (Fig. 4). It is possible 
that the more rapid growth of S in this case arises from 
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